import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

m = [1.33e20,3.98e14,4.9e12]
x = np.array([0,1.5e11,1.5e11+3.8e8])
y = np.array([0.0,0,0])
u = np.array([0.0,0,0])
v = np.array([0,2.88e4,2.88e4+1.02e3])

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(xlim=(-2e11,2e11),ylim=(-2e11,2e11))
ax.grid()

trace0, = ax.plot([],[],'-', lw=0.5)
trace1, = ax.plot([],[],'-', lw=0.5)
trace2, = ax.plot([],[],'-', lw=0.5)

pt0, = ax.plot([x[0]],[y[0]] ,marker='o')
pt1, = ax.plot([x[0]],[y[0]] ,marker='o')
pt2, = ax.plot([x[0]],[y[0]] ,marker='o')

k_text = ax.text(0.05,0.85,'',transform=ax.transAxes)
textTemplate = 't = %.3f days\n'

N = 1000
dt = 36000
ts =  np.arange(0,N*dt,dt)/3600/24
xs,ys = [],[]
for _ in ts:
    x_ij = (x-x.reshape(3,1))
    y_ij = (y-y.reshape(3,1))
    r_ij = np.sqrt(x_ij**2+y_ij**2)
    for i in range(3):
        for j in range(3):
            if i!=j :
                u[i] += (m[j]*x_ij[i,j]*dt/r_ij[i,j]**3)
                v[i] += (m[j]*y_ij[i,j]*dt/r_ij[i,j]**3)
    x += u*dt
    y += v*dt
    xs.append(x.tolist())
    ys.append(y.tolist())

xs = np.array(xs)
ys = np.array(ys)

def animate(n):
    trace0.set_data(xs[:n,0],ys[:n,0])
    trace1.set_data(xs[:n,1],ys[:n,1])
    #绘图时的地月距离扩大100倍，否则看不清
    tempX2S = xs[:n,1]+100*(xs[:n,2]-xs[:n,1])
    tempY2S = ys[:n,1]+100*(ys[:n,2]-ys[:n,1])
    trace2.set_data(tempX2S,tempY2S)

    pt0.set_data([xs[n,0]],[ys[n,0]])
    pt1.set_data([xs[n,1]],[ys[n,1]])
    tempX = xs[n,1]+100*(xs[n,2]-xs[n,1])
    tempY = ys[n,1]+100*(ys[n,2]-ys[n,1])
    pt2.set_data([tempX],[tempY])

    k_text.set_text(textTemplate % ts[n])
    return trace0, trace1, trace2, pt0, pt1, pt2, k_text

ani = animation.FuncAnimation(fig, animate, 
    range(N), interval=10, blit=True)

plt.show()
ani.save("3.gif")
